/* -------------------------------------------------------------------------- *
 *                      SimTK Core: SimTK Molmodel                            *
 * -------------------------------------------------------------------------- *
 * This is part of the SimTK Core biosimulation toolkit originating from      *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
 * Portions copyright (c) 2006-7 Stanford University and the Authors.         *
 * Authors: Samuel Flores                                                     *
 * Contributors:                                                              *
 *                                                                            *
 * Permission is hereby granted, free of charge, to any person obtaining a    *
 * copy of this software and associated documentation files (the "Software"), *
 * to deal in the Software without restriction, including without limitation  *
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
 * and/or sell copies of the Software, and to permit persons to whom the      *
 * Software is furnished to do so, subject to the following conditions:       *
 *                                                                            *
 * The above copyright notice and this permission notice shall be included in *
 * all copies or substantial portions of the Software.                        *
 *                                                                            *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
 * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
 * -------------------------------------------------------------------------- */

using namespace SimTK;
using namespace std;
const int maxChiWelds = 1000;
const int maxChains = 11;
const int maxChars    = 3000;
const int maxResidues = maxChars;
//#include "DistanceMatrix.h"
//#include <direct.h> // for getcwd
#include <stdlib.h>// for MAX_PATH
// #include "/Users/samuelflores/svn/simbody/src/ForceImpl.h" 
#include "SimTKmolmodel.h"
#include "Simbody.h"
#include "PeriodicScrubber.h"  
//#include "HarmonicVelocityRescalingThermostat.h"

//#include <VTKEventReporter.h>
#include <iostream>
#include <vector>

//#include "molmodel/internal/PeriodicVmdReporter.h"
#include "PeriodicPdbAndEnergyWriter.h"
#include "PeriodicPdbAndEnergySingleFrameWriter.h"
//#include "RNAToolbox/include/PeriodicPdbAndEnergyWriter.h"              
//#include "Water.h"    
#include "BaseInteractionParameterReader.h"
#include "ParameterReader.h"
#include "WaterDroplet.h"    
#include "TwoTransformForces.h"
#include "PeriodicParameterReaderUpdater.h"
//#include "BasePairingBonds.h"              
#include "SetChiMobility.h"              
#include "RigidifyFormedHelices.h"
//#include "ParameterReader2.h"
//#include "StackingForceSubsystem.h"

#include "AddHuntCrossleySpheres.h"
#include "AddNASTForces.h"
#include "AddBackboneOxygenForces.h"
#include "RandomizeRNACoordinates.h"
#include "KBBackboneTorsionForce.h"
#include "WadleyKeatingDuartePyleTorsionForce.h"
#include "SetSingleBondMobility.h" 
//using std::cout;
//using std::endl;


#include "release_joint.h"
//#include "lock_joint.h"





float PointToPlaneDistance (Vec3 Point1, Vec3 Normal1, Vec3 Point2) {
   float d = - dot(Point1,Normal1);
   float D = (d + dot(Point2,Normal1))/Normal1.norm();
   return D;
};
/*
class RepulsiveForce : public Force::Custom::Implementation { 
public: 
    RepulsiveForce(SimbodyMatterSubsystem& matter) : matter(matter) { 
    } 
    void calcForce(const State& state, Vector_<SpatialVec>& bodyForces,  
            Vector_<Vec3>& particleForces, Vector& mobilityForces) const { 
        for (MobilizedBodyIndex i(0); i < matter.getNBodies(); i++) { 
            const MobilizedBody& body1 = matter.getMobilizedBody(i); 
            for (MobilizedBodyIndex j(0); j < i; j++) { 
                const MobilizedBody& body2 = matter.getMobilizedBody(j); 
                Vec3 r = body1.getBodyOriginLocation(state)- 
                         body2.getBodyOriginLocation(state); 
                Real distance = r.norm(); 
                Vec3 force = 100*r/(distance*distance*distance); // 1 was 100, reduced for second tertiary contact formation
                bodyForces[i][1] += force; 
                bodyForces[j][1] -= force; 
            } 
        } 
    } 
    Real calcPotentialEnergy(const State& state) const { 
        double energy = 0.0; 
        for (MobilizedBodyIndex i(0); i < matter.getNBodies(); i++) { 
            const MobilizedBody& body1 = matter.getMobilizedBody(i); 
            for (MobilizedBodyIndex j(0); j < i; j++) { 
                const MobilizedBody& body2 = matter.getMobilizedBody(j); 
                Vec3 r = body1.getBodyOriginLocation(state)- 
                         body2.getBodyOriginLocation(state); 
                energy -= 1.0/r.norm(); 
            } 
        } 
        return energy; 
    } 
    bool dependsOnlyOnPositions() const { 
        return true; 
    } 
private: 
    SimbodyMatterSubsystem& matter; 
}; 
*/


/**
 * /brief This utility is used to get the name of an atom involved in hydrogen bonding. It was written to help with the problem of pulling together two halves of an A-form double helix. 
 *
 * By convention the first atom of the first residue is a bond donor "H", the second atom of the first residue is an accepter "O".
 * the opposite convention holds for the second residue.
 *
 * /param PdbResidueName is the type of base (AUCG).  Based on this parameter we will determine the name of the appropriate hydrogen-bonding atom.
 * 
 */

// this polymorphism soon to be obsolete

string ReturnHBondingAtomName(String myPdbResidueName, int firstOrSecondResidue, int firstOrSecondBond ) {
//
//
    assert((firstOrSecondResidue == 1 ) || (firstOrSecondResidue == 2));
    assert((firstOrSecondBond    == 1 ) || (firstOrSecondBond    == 2));
    string myG("G  ");
    string myC("C  ");
    string myA("A  ");
    string myU("U  ");

     cout<<"about to compare actual residue name = >"<<myPdbResidueName<< "< to "<<myG<<endl;
    if (!myG.compare(myPdbResidueName)) {
	 cout<<"This residue is a G"<<endl;
    	if ((firstOrSecondResidue == 1) ^ (firstOrSecondBond == 2)) {return "/H22";} else {return "/O6";}
	} else { cout<<"First residue is not a G"<<endl;} 
    if (!myC.compare(myPdbResidueName)) {
	 cout<<"This  residue is a C"<<endl;
    	if ((firstOrSecondResidue == 1) ^ (firstOrSecondBond == 2)) {return "/H41";} else {return "/O2";}
	} else { cout<<"First residue is not a C"<<endl;} 
    if (!myA.compare(myPdbResidueName)) {
	 cout<<"This  residue is a A"<<endl;
    	if ((firstOrSecondResidue == 1) ^ (firstOrSecondBond == 2)) {return "/H62";} else {return "/N1";}
	} else { cout<<"First residue is not a A"<<endl;} 
    if (!myU.compare(myPdbResidueName)) {
	 cout<<"This  residue is a U"<<endl;
    	if ((firstOrSecondResidue == 1) ^ (firstOrSecondBond == 2)) {return "/H3" ;} else {return "/O2";}
	} else { cout<<"First residue is not a U"<<endl;} 

}
/**
 * 
 * 
 * /param 
 * myPdbResidueName1,2 must be one of "A","C","G","U".
 * bondingEdge1,2 must be one of "WatsonCrick","Hoogsteen","Sugar","Bifurcated".
 * glycosidicBondOrientation must be either "Cis" or "Trans".
 *
 */

/**
 * /brief This utility is used for the Concentricity constraint.  In this constraint, a certain atom on residue 1 is connected by a LinearSpring to a certain point OUTSIDE of an atom on residue 2.  location2 is that point, in the body frame of the atom on residue 2.
 * 
 */

class  SetPolynucleotideChiBondMobility  { public:

    SetPolynucleotideChiBondMobility (
        Biopolymer & myChain,
	int startResidueNumber,
	int endResidueNumber,
	LeontisWesthofClass  & myLeontisWesthofClass ,    
        SimbodyMatterSubsystem & matter,
        State& state,
        Constraint myWeld1[maxChiWelds],
        Constraint myWeld2[maxChiWelds],
        Constraint myWeld3[maxChiWelds]

	)   
	{   
        for (int q = (startResidueNumber); q <= (endResidueNumber); q++) {
            SetChiMobility(myChain,q,myLeontisWesthofClass , matter,state,myWeld1[q],myWeld2[q],myWeld3[q]    );
            //SetChiMobility(myChain,q,myLeontisWesthofClass , mobility    );
        }
        }
    SetPolynucleotideChiBondMobility (
        Biopolymer & myChain,
	int startResidueNumber,
	int endResidueNumber,
	LeontisWesthofClass  & myLeontisWesthofClass ,    
	BondMobility::Mobility mobility = BondMobility::Free   
	)   
	{   
        for (int q = (startResidueNumber); q <= (endResidueNumber); q++) {
            //SetChiMobility(myChain,q,myLeontisWesthofClass , matter,state    );
             //cout<<"[Repel.h:SetPolynucleotideChiBondMobility] setting residue :"<<q<<" to BondMobility :"<<mobility<<endl;
            SetChiMobility(myChain,q,myLeontisWesthofClass , mobility    );
        }
        }
};       
class CentrifugalForce : public Force::Custom::Implementation { 
public: 
    CentrifugalForce(SimbodyMatterSubsystem& matter) : matter(matter) {  
    } 
    void calcForce(const State& state, Vector_<SpatialVec>& bodyForces,  
            Vector_<Vec3>& particleForces, Vector& mobilityForces) const { 
        for (MobilizedBodyIndex i(0); i < matter.getNBodies(); i++) { 
            const MobilizedBody& body1 = matter.getMobilizedBody(i); 
            Vec3 r = body1.getBodyOriginLocation(state); 
            Vec3 force = r;//(distance*distance*distance); 
	    // cout<<"r,force ="<<r<<","<<force<<endl;
            bodyForces[i][1] += force; 
         } 
    } 
    Real calcPotentialEnergy(const State& state) const { 
        double energy = 0.0; 
        for (MobilizedBodyIndex i(0); i < matter.getNBodies(); i++) { 
            const MobilizedBody& body1 = matter.getMobilizedBody(i); 
            Vec3 r = body1.getBodyOriginLocation(state); 
            energy += (r.norm())*(r.norm());//(distance*distance*distance); 
        } 
	 cout<<"just computed energy of "<<energy<<endl; 
        return energy; 
    } 
    bool dependsOnlyOnPositions() const { 
        return true; 
    } 
private: 
    SimbodyMatterSubsystem& matter; 
}; 
    
 

class    ConstrainedDynamics : public Compound {
public:
    LeontisWesthofClass myLeontisWesthofClass;
    int setDefaults() {
    }
    Compound runDynamics(ParameterReader & myParameterReader) {
	// cout <<" myParameterReader.helixBondMobility ="<<myParameterReader.helixBondMobility<<endl;
        myLeontisWesthofClass.initialize((myParameterReader.leontisWesthofInFileName));
        CompoundSystem system;
        SimbodyMatterSubsystem  matter(system);
        GeneralForceSubsystem forces(system);
        GeneralContactSubsystem contacts(system);

        DuMMForceFieldSubsystem backboneOxygenForces;
	    AddBackboneOxygenForces myAddBackboneOxygenForces;
        myAddBackboneOxygenForces.initialize(backboneOxygenForces,myParameterReader);//(myMagnesiumIonVec[0])*,myParameterReader); 
        DuMMForceFieldSubsystem nastForces;
        DuMMForceFieldSubsystem wkdpForces;
	if (myParameterReader.addNASTForces) {
	    system.addForceSubsystem( static_cast<SimTK::ForceSubsystem&> (nastForces) );
        }
        if (myParameterReader.wkdpGlobalBondTorsionScaleFactor > 0.0) {
            system.addForceSubsystem( static_cast<SimTK::ForceSubsystem&> (wkdpForces) );
        }
	//if (myParameterReader.addBackboneOxygenForces) {
	    system.addForceSubsystem( static_cast<SimTK::ForceSubsystem&> (backboneOxygenForces) );
            //VanderWallSphere myVanderWallSphere( forces,  backboneOxygenForces,Vec3(0),myParameterReader.vanderWallSphereRadius,.1,1.0); 
            //OrientedBoundingBox myOrientedBoundingBox(Transform(Rotation(0,Vec3(0,0,1)),Vec3(0)), Vec3(100,100,100));
	    
        //}
	//DuMMForceFieldSubsystem nastForces(system);
        DuMMForceFieldSubsystem dumm(system);
        //DuMMForceFieldSubsystem dumm(system);
        HuntCrossleyContact myHuntCrossleyContact(system);
	//scf
	//nastForces.loadAmber99Parameters ();

        ContactSetIndex contactSet = contacts.createContactSet();
        ContactSetIndex contactSetLargeSpheres = contacts.createContactSet();
        HuntCrossleyForce hc(forces, contacts, contactSet);
        HuntCrossleyForce hcLargeSpheres(forces, contacts, contactSetLargeSpheres);

        filebuf fb;
        fb.open ((myParameterReader.tinkerParameterFileName).c_str(),ios::in);
        istream is(&fb);
	if (myParameterReader.loadTinkerParameterFile) {
            cout<<"You have specified tinkerParameterFileName = "<<myParameterReader.tinkerParameterFileName<<".  Checking this file.."<<endl;
            SimTK_ERRCHK_ALWAYS(fb.is_open(),"[Repel.h]", "The Tinker parameter file you specified could not be opened.  Please check your tinkerParameterFileName parameter, or set \"loadTinkerParameterFile 0\" to use the hard-coded Tinker parameters instead.");
            dumm.populateFromTinkerParameterFile (is);}
        else 
            dumm.loadAmber99Parameters();
        fb.close();
        if (myParameterReader.kbBackboneTorsionGlobalScaleFactor > 0.0) {
            DefineKBBackboneTorsion(dumm,myParameterReader);
        }
        if (myParameterReader.wkdpGlobalBondTorsionScaleFactor > 0.0) {
            DefineWKDPTorsion(wkdpForces,myParameterReader);
        }
        //dumm.populateFromTinkerParameterFile (std::istream("/Users/samuelflores/svn/RNAToolbox/include/resources/tinker_amber99_sam.prm"));
        dumm.setCoulombGlobalScaleFactor(myParameterReader.globalCoulombScaleFactor);
        dumm.setBondTorsionGlobalScaleFactor(myParameterReader.globalBondTorsionScaleFactor);
        dumm.setGbsaGlobalScaleFactor(myParameterReader.globalGbsaScaleFactor);
        dumm.setVdwGlobalScaleFactor(myParameterReader.globalVdwScaleFactor);
        dumm.setBondStretchGlobalScaleFactor(myParameterReader.globalBondStretchScaleFactor);
        dumm.setBondBendGlobalScaleFactor(myParameterReader.globalBondBendScaleFactor);
        dumm.setAmberImproperTorsionGlobalScaleFactor(myParameterReader.globalAmberImproperTorsionScaleFactor);
        if (myParameterReader.kbBackboneTorsionGlobalScaleFactor == 0) dumm.setCustomBondTorsionGlobalScaleFactor(0); else dumm.setCustomBondTorsionGlobalScaleFactor(1);
        dumm.setUseMultithreadedComputation(myParameterReader.useMultithreadedComputation);
        //dumm.setAllGlobalScaleFactors(0);

	map<const string,String>::iterator sequencesIterator = (myParameterReader.sequences).begin();
        Biopolymer myMolecule[maxChains];
        //vector<Biopolymer> myMolecule(maxChains);
        Biopolymer myProtein[maxChains];
        for (int i = 0; i < (myParameterReader.sequences).size(); i++) 
            myMolecule[i] = RNA(myParameterReader.sequences[(myParameterReader.chainId[i]).c_str()] ,1);
        for (int i = 0; i < (myParameterReader.proteinSequences).size(); i++) 
            myProtein[i] = Protein(myParameterReader.proteinSequences[(myParameterReader.proteinChainId[i]).c_str()] );
	//if (myParameterReader.verbose) cout<<"sequence set to ="<<myParameterReader.sequences[(*sequencesIterator).first]<<endl;
        vector<MagnesiumIon> myMagnesiumIonVec(myParameterReader.numMagnesiumIons);            

        time_t rawtime;
        struct tm * timeinfo;
        time ( &rawtime );
        timeinfo = localtime ( &rawtime );
        if (myParameterReader.verbose) cout<<"[Repel.h] check -2: time:"<< asctime (timeinfo)   <<endl;

        //if (setRepulsiveForce) Force::Custom(forces, new RepulsiveForce(matter)); 
        //assert (!((myParameterReader.readPreviousFrameFile == 1) && (readInQVector == 1))); 
	//map<const char*, int   , strCmp>::iterator firstResidueNumbersIterator = myParameterReader.firstResidueNumbersIterator;
        for (int i = 0; i < (myParameterReader.sequences).size(); i++) {myMolecule[i].renumberPdbResidues( myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()) );}
        for (int i = 0; i < (myParameterReader.proteinSequences).size(); i++) {myProtein[i].renumberPdbResidues( myParameterReader.getProteinFirstResidueNumbers((myParameterReader.proteinChainId[i]).c_str()) -1);} // -1 accounts for the ACE residue chris adds at the beginning of the chain
	//cout<<"Renumbering starting at :"<<myParameterReader.firstResidueNumbers[(*firstResidueNumbersIterator).first]<<endl;
//	for (int q = 0; q<myParameterReader.NResidues;q++) for (int i = 0; i < (myParameterReader.sequences).size(); i++) (myMolecule[i].updResidue(Compound::Index(q))).setPdbResidueNumber(myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str())+q);
        for (int i = 0; i < (myParameterReader.numMagnesiumIons); i++) {
	    myMagnesiumIonVec[i].setPdbChainId(*((myParameterReader.magnesiumIonChainId).c_str()));
	    myMagnesiumIonVec[i].setPdbResidueNumber(i);
        }
        for (int i = 0; i < (myParameterReader.sequences).size(); i++) 
            {
	    myMolecule[i].setPdbChainId(*(((myParameterReader.chainId[i]).substr(0,1)).c_str()));
	    }
        for (int i = 0; i < (myParameterReader.proteinSequences).size(); i++) 
            {
	    myProtein[i].setPdbChainId(*(((myParameterReader.proteinChainId[i]).substr(0,1)).c_str()));
	    }
        if (myParameterReader.readPreviousFrameFile) 
            for (int i = 0; i < (myParameterReader.sequences).size(); i++) system.adoptCompound(myMolecule[i],Vec3(  0,0,0));
        else
            for (int i = 0; i < (myParameterReader.sequences).size(); i++) system.adoptCompound(myMolecule[i],Vec3(i*2 ,0,0));
        for (int i = 0; i < (myParameterReader.proteinSequences).size(); i++) system.adoptCompound(myProtein[i],Vec3(i*4 ,0,0));
        //float myParameterReader.magnesiumIonRadius =25.0;
        float magnesiumIonHelixAdvance =10.0;
        float magnesiumIonAngleIncrement = 20.0;
        for (int i = 0; i < (myParameterReader.numMagnesiumIons); i++) {
            system.adoptCompound(myMagnesiumIonVec[i],Vec3(myParameterReader.magnesiumIonRadius*cos(i*magnesiumIonAngleIncrement/Rad2Deg),myParameterReader.magnesiumIonRadius*sin(i*magnesiumIonAngleIncrement/Rad2Deg),i*magnesiumIonAngleIncrement/360*magnesiumIonHelixAdvance ));
        if (myParameterReader.waterDropletMake) 
            WaterDroplet myWaterDroplet(system, dumm, forces,matter, myParameterReader,contacts, contactSet);
        }




        if (myParameterReader.readPreviousFrameFile == 1) {

            cout<<"Fitting structure to "<<myParameterReader.previousFrameFileName<<endl;
            ifstream previousFrameFile (myParameterReader.previousFrameFileName.c_str());
            assert(previousFrameFile.good());
            PdbStructure pdbStructure(previousFrameFile);
            if (myParameterReader.verbose) pdbStructure.write(std::cout);
            if (myParameterReader.readMagnesiumPositionsFromFile) for (int i = 0; i < (myParameterReader.numMagnesiumIons); i++) {
                Compound::AtomTargetLocations magnesiumAtomTargets = myMagnesiumIonVec[i].createAtomTargets(pdbStructure); 
                if (myParameterReader.verbose) cout<<"[Repel.h] magensiumAtomTargets: "<<magnesiumAtomTargets[AtomIndex(0)]<<endl;
                if (myParameterReader.verbose) cout<<"[Repel.h] magensiumAtomTargets: "<<magnesiumAtomTargets.size()       <<endl;
		myMagnesiumIonVec[i].fitDefaultConfiguration(magnesiumAtomTargets,myParameterReader.fitDefaultTolerance);
                //myMagnesiumIonVec[i].matchDefaultTopLevelTransform(magnesiumAtomTargets);
            }
            //assert(0);

        for (int i = 0; i < (myParameterReader.sequences).size(); i++) {
             
           
            Compound::AtomTargetLocations atomTargets = myMolecule[i].createAtomTargets(pdbStructure); 
	    if (myParameterReader.verbose) cout<<"size of atomTargets ="<<atomTargets.size()<<endl;	
	    if (myParameterReader.verbose) cout<<"number of atoms in myMolecule[0] ="<<myMolecule[0].getNAtoms()<<endl;

            {
            //myMolecule[i].matchDefaultConfiguration(atomTargets ); //planarity tolerance is in Radians, if Sherm's email is to be believed
            time ( &rawtime );
            timeinfo = localtime ( &rawtime );
            if (myParameterReader.verbose) cout<<"[Repel.h] check -1a: time:"<< asctime (timeinfo)   <<endl;
            myMolecule[i].matchDefaultAtomChirality(atomTargets,.010);
            time ( &rawtime );
            timeinfo = localtime ( &rawtime );
            if (myParameterReader.verbose) cout<<"[Repel.h] check -1b: time:"<< asctime (timeinfo)   <<endl;
            myMolecule[i].matchDefaultBondLengths(atomTargets);
            time ( &rawtime );
            timeinfo = localtime ( &rawtime );
            if (myParameterReader.verbose) cout<<"[Repel.h] check -1c: time:"<< asctime (timeinfo)   <<endl;
            myMolecule[i].matchDefaultBondAngles(atomTargets);
            time ( &rawtime );
            timeinfo = localtime ( &rawtime );
            if (myParameterReader.verbose) cout<<"[Repel.h] check -1d: time:"<< asctime (timeinfo)   <<endl;
            myMolecule[i].matchDefaultDihedralAngles(atomTargets);
            time ( &rawtime );
            timeinfo = localtime ( &rawtime );
            if (myParameterReader.verbose) cout<<"[Repel.h] check -1e: time:"<< asctime (timeinfo)   <<endl;
            if (! myParameterReader.matchDefaultSkipTopLevelTransform) myMolecule[i].matchDefaultTopLevelTransform(atomTargets);
            time ( &rawtime );
            timeinfo = localtime ( &rawtime );
            if (myParameterReader.verbose) cout<<"[Repel.h] check -1f: time:"<< asctime (timeinfo)   <<endl;
            cout<<"[Repel.h] About to fit RNA chain "<<myParameterReader.chainId[i]<<".  "<<endl;//If this step fails your fitDefaultTolerance parameter is probably too small."<<endl;
            }
            /*
            if (atomTargets.size() > 0)  {
                cout<<"[Repel.h] atomTargets.size() ="<<atomTargets.size()<<endl;
		for (int k = 0 ; k < atomTargets.size() ; k++) cout<<"[Repel.h] k, atomTargets[k] :"<<k<<","<<atomTargets[Compound::AtomIndex(k)]<<","<<myMolecule[i].getAtomName(Compound::AtomIndex(k))<<endl;
                cout<<"[Repel.h] myMolecule[i].getNAtoms = "<<(myMolecule[i]).getNAtoms()<<endl;

                //myMolecule[i].setCompoundBondMobility(BondMobility::Free   );

   
                    if (0) if (myParameterReader.setChiBondMobility) {
                        SetPolynucleotideChiBondMobility (
                            myMolecule[i],           
	                    0,
        		    myMolecule[i].getNResidues()-1,//(myParameterReader.sequences[i]).length()-1,// - myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str())),
	                    //(myParameterReader.secondBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str())),
	                    myLeontisWesthofClass,
	                    myParameterReader.chiBondMobility   
	                );
                 
                    }
 	        myMolecule[i].fitDefaultConfiguration(atomTargets,myParameterReader.fitDefaultTolerance);
		if (myParameterReader.verbose) cout<<"[Repel.h] Just executed matchDefault and fitDefaultConfiguration."<<endl;
            }
            */
            //myMolecule[i].matchDefaultConfiguration(atomTargets,Compound::Match_Idealized); //planarity tolerance is in Radians, if Sherm's email is to be believed
								       //planarityTolerance should definitely be a large number like 1.57.  otherwise base bonding will get out of whack.
	    previousFrameFile.close();
	    	
        }
        time ( &rawtime );
        timeinfo = localtime ( &rawtime );
        if (myParameterReader.verbose) cout<<"[Repel.h] check 0: time:"<< asctime (timeinfo)   <<endl;

        for (int i = 0; i < (myParameterReader.proteinSequences).size(); i++) {
            if (myParameterReader.verbose) cout<<"Fitting structure to "<<myParameterReader.previousFrameFileName<<endl;
            ifstream previousFrameFile (myParameterReader.previousFrameFileName.c_str());
	    cout<<"Checking for "<<myParameterReader.previousFrameFileName<<endl;
            assert(previousFrameFile.good());
            PdbStructure pdbStructure(previousFrameFile);
            Compound::AtomTargetLocations atomTargets = myProtein[i].createAtomTargets(pdbStructure); 
	    if (myParameterReader.verbose) cout<<"size of atomTargets ="<<atomTargets.size()<<endl;	
	    cout<<"number of atoms in myMolecule[0] ="<<myMolecule[0].getNAtoms()<<endl;
            myProtein[i].matchDefaultConfiguration(atomTargets); //planarity tolerance is in Radians, if Sherm's email is to be believed
								       //planarityTolerance should definitely be a large number like 1.57.  otherwise base bonding will get out of whack.
	    previousFrameFile.close();
            myProtein[i].writeDefaultPdb("afterMatchDefault.pdb",Vec3(0));
	    	
        }
	}
		


	if (myParameterReader.setOverallBondMobility) for (int i = 0; i < (myParameterReader.sequences).size(); i++) {
            myMolecule[i].setCompoundBondMobility(myParameterReader.overallBondMobility);
        }
        time ( &rawtime );
        timeinfo = localtime ( &rawtime );
        if (myParameterReader.verbose) cout<<"[Repel.h] check -2a: time:"<< asctime (timeinfo)   <<endl;
        if (myParameterReader.piecewiseRigidify) { 
	    //cout<<"[Repel.h] : myParameterReader.piecewiseRigidify .. myParameterReader.numRigidSegments[(myParameterReader.chainId[0]).c_str()][0], myParameterReader.rigidSegmentStartNumber[0], myParameterReader.rigidSegmentEndNumber[0] ="<<  myParameterReader.numRigidSegments[(myParameterReader.chainId[0]).c_str()]<<","<< myParameterReader.rigidSegmentStartNumber[0]<<","<< myParameterReader.rigidSegmentEndNumber[0]<<endl; 
        for (int i = 0; i < (myParameterReader.sequences).size(); i++) for (int q=0; q < (myParameterReader.numRigidSegments[(myParameterReader.chainId[i]).c_str()]); q++){
	    if (myParameterReader.verbose) cout<<"[Repel.h] : myParameterReader.piecewiseRigidify .. myParameterReader.numRigidSegments[(myParameterReader.chainId[i]).c_str()] , myParameterReader.rigidSegmentStartNumber[q][i], myParameterReader.rigidSegmentEndNumber[q][i] ="<<  myParameterReader.numRigidSegments[(myParameterReader.chainId[i]).c_str()]<<","<< myParameterReader.rigidSegmentStartNumber[q][i]<<","<< myParameterReader.rigidSegmentEndNumber[q][i]<<endl; 
            RNA& rna = static_cast<RNA&>(myMolecule[i]); 
            if (myParameterReader.verbose) cout<<"[Repel.h] myParameterReader.chainId[i]="<<myParameterReader.chainId[i]<<endl;
            if (myParameterReader.BPEdge[(myParameterReader.rigidSegmentStartNumber[q][i]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()))][i].compare("WatsonCrick") ==0 ){
                for (int j = 0; j < myParameterReader.numChains; j++) assert(myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[j]).c_str()) > 0);
                //cout<<"[Repel.h]  q (segment#) ,i (chainId) ,myParameterReader.helixBondMobility,(myParameterReader.rigidSegmentStartNumber[q][i]-myParameterReader.//firstResidueNumbers[(myParameterReader.chai    nId[i]).c_str()]),(myParameterReader.rigidSegmentEndNumber[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()))"<< q<<","<<i<<","<<myParameterReader.helixBondMobility<<","<<myParameterReader.rigidSegmentStartNumber[q]<<"-"<<myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str())<<","<<myParameterReader.rigidSegmentEndNumber[q]<<"-"<<myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str())<<endl;
		rna.setRNABondMobility(    myParameterReader.helixBondMobility,(myParameterReader.rigidSegmentStartNumber[q][i]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str())),(myParameterReader.rigidSegmentEndNumber[q][i]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()))); 
                }
            else {
		//set myParameterReader.loopBondMobility.  cut this out for debugging, might want it back later
		if (myParameterReader.setLoopBondMobility) (rna).setRNABondMobility(myParameterReader.loopBondMobility,(myParameterReader.rigidSegmentStartNumber[q][i]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str())),(myParameterReader.rigidSegmentEndNumber[q][i]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()))); 
                }
            myMolecule[i] = rna;

        }

	}
 

        for (int i = 0; i < (myParameterReader.sequences).size(); i++)  for (int q=0;q<myParameterReader.numRigidSegments[(myParameterReader.chainId[i]).c_str()]; q++)         {


            if (myParameterReader.setChiBondMobility) {
		if (myParameterReader.verbose) cout<<"[Repel.h] myParameterReader.setChiBondMobility ="<<myParameterReader.setChiBondMobility<<endl;	
                //for (int i = 0; i < (myParameterReader.sequences).size(); i++) 
		for (int k = 0; k < myParameterReader.numChains; k++) assert(myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[k]).c_str()) > 0);
                if (myParameterReader.verbose) cout<<"[Repel.h]  myParameterReader.rigidSegmentStartNumber[q][i] , myParameterReader.rigidSegmentEndNumber[q][i] "<<myParameterReader.rigidSegmentStartNumber[q][i] <<","<< myParameterReader.rigidSegmentEndNumber[q][i]<<   endl;
                SetPolynucleotideChiBondMobility (
                    myMolecule[i],           
	            (myParameterReader.rigidSegmentStartNumber[q][i]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str())),
	            (myParameterReader.rigidSegmentEndNumber[q][i]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str())),
	            myLeontisWesthofClass,
	            myParameterReader.chiBondMobility   
	            );
                }
 
        //Constraint residueWeld1[maxChiWelds];//sized according to the maximum possible number of Welds we will need.
        }

        time ( &rawtime );
        timeinfo = localtime ( &rawtime );
        if (myParameterReader.verbose) cout<<"[Repel.h] check -2b: time:"<< asctime (timeinfo)   <<endl;
//scf added
        //vector<Biopolymer> myMoleculeVector;//(myMolecule, myMolecule  + sizeof(myMolecule ) / sizeof(Biopolymer) );
	//for (int i = 0; i<maxChains; i++) myMoleculeVector[i] = myMolecule[i];
	for (int q=0;q<myParameterReader.numBasePairs;q++) {
            int i = -11111;
            int j = -11111;
            for (int k = 0; k < myParameterReader.numChains; k++) {
                if (strcmp((myParameterReader.chainId[k]).c_str(),(myParameterReader.firstBPChain[q]).c_str()) == 0) i = k;
                if (strcmp((myParameterReader.chainId[k]).c_str(),(myParameterReader.secondBPChain[q]).c_str()) == 0) j = k;
            }   
            assert(i>=0);
            assert(j>=0);
	    //assert(i==j);
            /*SimTK_ERRCHK_ALWAYS(
                    ((i == j)) ,
                    "[Repel.h]",
                    "For the basePair Rigid and the SingleBondMobility parameters you must have the start and end points on the same chain."); */
        if ((myParameterReader.firstBPEdge[q]).compare("Rigid") == 0) 
           {

                if ((myParameterReader.firstBPEdge[q]).compare("Rigid") == 0) {
                    if (myParameterReader.verbose) cout<<"[Repel.h])rameterReader.firstBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()) , myParameterReader.secondBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[j]).c_str()) :"<< myParameterReader.firstBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()) <<","<< myParameterReader.secondBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[j]).c_str())  <<endl; 
                    if (myParameterReader.verbose) cout<<"[Repel.h]myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()) , myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()) :"<<myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()) <<","<< myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[j]).c_str())  <<endl; 
                   (static_cast<RNA&>(myMolecule[i])).setRNABondMobility(BondMobility::Rigid ,myParameterReader.firstBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()) , myParameterReader.secondBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str())  );
   
                    if (myParameterReader.setChiBondMobility) {
                        assert (i==j);
                        SetPolynucleotideChiBondMobility (
                            myMolecule[i],           
	                    (myParameterReader.firstBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str())),
	                    (myParameterReader.secondBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str())),
	                    myLeontisWesthofClass,
	                    myParameterReader.chiBondMobility   
	                );
                 
                    }


		}

           } 
           else if ((myParameterReader.firstBPEdge[q]).compare("SingleBondMobilityRigid") == 0)
               SetSingleBondMobility(myMolecule[i],
                   myParameterReader.firstBPResidue[q]- myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()),
                   myParameterReader.secondBPResidue[q]- myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()),
                   myParameterReader.secondBPEdge[q],
                   myParameterReader.orientationBP[q],
                   BondMobility::Rigid);
           else if ((myParameterReader.firstBPEdge[q]).compare("SingleBondMobilityTorsion")== 0)
               SetSingleBondMobility(myMolecule[i],
                   myParameterReader.firstBPResidue[q]- myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()),
                   myParameterReader.secondBPResidue[q]- myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()),
                   myParameterReader.secondBPEdge[q],
                   myParameterReader.orientationBP[q],
                   BondMobility::Torsion);
           else if ((myParameterReader.firstBPEdge[q]).compare("SingleBondMobilityFree" ) == 0)
               SetSingleBondMobility(myMolecule[i],
                   myParameterReader.firstBPResidue[q]- myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()),
                   myParameterReader.secondBPResidue[q]- myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()),
                   myParameterReader.secondBPEdge[q],
                   myParameterReader.orientationBP[q],
                   BondMobility::Free );
           }


        int t  = 0;
        //for (int q=0;q<myParameterReader.numRigidSegments; q++) 

		//(static_cast<RNA&>(myMolecule[0])).setRNABondMobility(BondMobility::Rigid,0,0);
  //       (static_cast<RNA&>(myMolecule[0])).setRNABondMobility(BondMobility::Rigid,1,1);
  //      (static_cast<RNA&>(myMolecule[0])).setRNABondMobility(BondMobility::Rigid,2,2);
		// (static_cast<RNA&>(myMolecule[0])).setRNABondMobility(BondMobility::Rigid,3,3);//Additionalline
		// (static_cast<RNA&>(myMolecule[0])).setRNABondMobility(BondMobility::Rigid,4,4);
		// (static_cast<RNA&>(myMolecule[0])).setRNABondMobility(BondMobility::Rigid,5,5);
		// (static_cast<RNA&>(myMolecule[0])).setRNABondMobility(BondMobility::Rigid,6,6);
		//
		
		 

        if (myParameterReader.rigidifyFormedHelices) {
        //for (int i = 0; i < (myParameterReader.sequences).size(); i++) 
	{
            //RNA& rna = static_cast<RNA&>(myMolecule[i]);
	    RigidifyFormedHelices(myParameterReader, myMolecule ,myParameterReader.helixBondMobility);
	    //myMolecule[i] = rna;	
        }}
        time ( &rawtime );
        timeinfo = localtime ( &rawtime );
        if (myParameterReader.verbose) cout<<"[Repel.h] check -1: time:"<< asctime (timeinfo)   <<endl;

        if (myParameterReader.verbose) {
        //const State & state5 = state;
        string myFileName5 = "check.4.pdb";
        ofstream  myOfstream5(  myFileName5.c_str());
	for (int i = 0; i < (myParameterReader.sequences).size(); i++) {
            myMolecule[i].writeDefaultPdb(( myOfstream5),Transform(Vec3(0)));
            if (myParameterReader.verbose) cout<<"[Repel.h] check 5 writing molecule i="<<i<<endl;   
        }
        time ( &rawtime );
        timeinfo = localtime ( &rawtime );
        if (myParameterReader.verbose) cout<<"[Repel.h] check 3: time:"<< asctime (timeinfo)   <<endl;
        }

        if (myParameterReader.rigidifyTermini) {
	    for (int i = 0; i < (myParameterReader.sequences).size(); i++) {
                RNA& rna = static_cast<RNA&>(myMolecule[i]);
                if (myParameterReader.verbose) cout <<"[Repel.h] myParameterReader.rigidifyTermini: myParameterReader.calcLowestBondingResidue(myParameterReader.chainId[i]), myParameterReader.calcHighestBondingResidue(myParameterReader.chainId[i]), myParameterReader.getFirstResidueNumbers(myParameterReader.chainId[i]), myParameterReader.sequences[myParameterReader.chainId[i]]).length() ="<< myParameterReader.calcLowestBondingResidue(myParameterReader.chainId[i])<<","<< myParameterReader.calcHighestBondingResidue(myParameterReader.chainId[i])<<","<< myParameterReader.getFirstResidueNumbers(myParameterReader.chainId[i])<<","<< (myParameterReader.sequences[myParameterReader.chainId[i]]).length()<<endl;
                rna.setRNABondMobility(    myParameterReader.helixBondMobility,0,myParameterReader.calcLowestBondingResidue(myParameterReader.chainId[i]) - myParameterReader.getFirstResidueNumbers(myParameterReader.chainId[i]) - 1 );
                rna.setRNABondMobility(    myParameterReader.helixBondMobility,myParameterReader.calcHighestBondingResidue(myParameterReader.chainId[i]) - myParameterReader.getFirstResidueNumbers(myParameterReader.chainId[i]) + 1 ,(myParameterReader.sequences[myParameterReader.chainId[i]]).length()-1 );
                myMolecule[i] = rna;
	    }	
	}

	for (int i = 0; i < (myParameterReader.sequences).size(); i++) myMolecule[i].writeDefaultPdb("after-matchdefault.pdb",Vec3(0));
        time ( &rawtime );
        timeinfo = localtime ( &rawtime );
        if (myParameterReader.verbose) cout<<"[Repel.h] check 1: time:"<< asctime (timeinfo)   <<endl;
        system.modelCompounds();
        time ( &rawtime );
        timeinfo = localtime ( &rawtime );
        if (myParameterReader.verbose) cout<<"[Repel.h] check 2: time:"<< asctime (timeinfo)   <<endl;
        //the following will not work if the C3* atoms are rigidified.  
        //turned this off when i found that with the debug Simbody libraries I was getting an invalid atomIx 0
        if (myParameterReader.addNASTForces) {
	    AddNASTForces myAddNASTForces;
            myAddNASTForces.initialize(nastForces,myParameterReader); 
            for (int i = 0; i < (myParameterReader.sequences).size(); i++) myAddNASTForces.addForces( system, matter , nastForces, myParameterReader, myMolecule[i]);
	}
        //scf
        if (myParameterReader.wkdpGlobalBondTorsionScaleFactor > 0) {
	    AddWKDPForces myAddWKDPForces;
            myAddWKDPForces.initialize(wkdpForces,myParameterReader); 
            for (int i = 0; i < (myParameterReader.sequences).size(); i++) myAddWKDPForces.addForces( system, matter , wkdpForces, myParameterReader, myMolecule[i]);
	}
        //if (myParameterReader.addBackboneOxygenForces)  // this is no longer a supported parameter.  Just turn Vdw and Coulomb on and off separately.
        //turned this off when i found that with the debug Simbody libraries I was getting an invalid atomIx 0
        if ((myParameterReader.backboneOxygenGlobalCoulombScaleFactor>0 ) || (myParameterReader.backboneOxygenGlobalVdwScaleFactor>0)){
            for (int i = 0; i < (myParameterReader.sequences).size(); i++) if (myMolecule[i].getNumResidues() > 1) myAddBackboneOxygenForces.addForces( system, matter , backboneOxygenForces, myParameterReader, myMolecule[i],1, myMolecule[i].getNumResidues()-1);
            for (int i = 0; i < (myParameterReader.numMagnesiumIons) ; i++) myAddBackboneOxygenForces.addForces( system, matter , backboneOxygenForces, myParameterReader, myMagnesiumIonVec[i]);
	}



        if (myParameterReader.verbose) {
        //const State & state5 = state;
        string myFileName5 = "check.4.5.pdb";
        ofstream  myOfstream5(  myFileName5.c_str());
	for (int i = 0; i < (myParameterReader.sequences).size(); i++) {
            myMolecule[i].writeDefaultPdb(( myOfstream5),Transform(Vec3(0)));
            if (myParameterReader.verbose) cout<<"[Repel.h] check 5 writing molecule i="<<i<<endl;   
        }
        time ( &rawtime );
        timeinfo = localtime ( &rawtime );
        if (myParameterReader.verbose) cout<<"[Repel.h] check 3: time:"<< asctime (timeinfo)   <<endl;
        }

        State &  state = system.updDefaultState();//;// = system.realizeTopology();
        //#system.modelCompounds();
        // This is the first Topology realization called.

		// This part is added to the code to define prescribed motions a*cos(bt) to the first residue.
		

     Vector_<Real> coefficients1(2); // two prescribed motions : coefficients1[0] * cos(coefficients1[1]*t)
	                                //                          coefficients1[2] * cos(coefficients1[3]*t)
    coefficients1[0] = 0.05; 
    coefficients1[1] = 20.0;
	//coefficients1[2] = 0.01;
	//coefficients1[3] = 30.0;
    Constraint::PrescribedMotion constraint0(matter,new Function::HarmonicCosShifted(coefficients1),MobilizedBodyIndex(1), MobilizerQIndex(4));
    
	//coefficients1[0] = 0.01; 
 //   coefficients1[1] = 0;
	//coefficients1[2] = 0;
	//
	//Constraint::PrescribedMotion constraint0(matter,new Function::Polynomial(coefficients1),MobilizedBodyIndex(1), MobilizerQIndex(4));
 //   
	
	Vector_<Real> coefficients2(2);
    coefficients2[0] = 0.05;
    coefficients2[1] = 40.0;
	//coefficients2[2] = 10.0;
	//coefficients2[3] = 15.0;
	//coefficients2[4] = 20.0;
 //   coefficients2[5] = 25.0;
	//coefficients2[6] = 30.0;
	//coefficients2[7] = 35.0;
	//coefficients2[8] = 40.0;
 //   coefficients2[9] = 45.0;
	//coefficients2[10] = 50.0;
	Constraint::PrescribedMotion constraint1(matter,new Function::HarmonicCosShifted(coefficients2),MobilizedBodyIndex(1), MobilizerQIndex(5));


	////Vector_<Real> coefficients3(4);
    //coefficients3[0] = 0.05;
    //coefficients3[1] = 20.0;
	//coefficients2[2] = 0.01;
	//coefficients2[3] = 30.0;
	//Constraint::PrescribedMotion constraint2(matter,new Function::HarmonicCosShifted(coefficients3),MobilizedBodyIndex(1), MobilizerQIndex(6));

	/*Vector_<Real> coefficients3(4);
    coefficients3[0] = 0.05;
    coefficients3[1] = 40.0;
	coefficients3[2] = 0.1;
	coefficients3[3] = 60.0;
	Constraint::PrescribedMotion constraint7(matter,new Function::HarmonicCosShifted(coefficients3),MobilizedBodyIndex(1), MobilizerQIndex(6));*/


	/*coefficients2[0] = 0.01; 
    coefficients2[1] = 0;
	coefficients2[2] = 0;
    Constraint::PrescribedMotion constraint1(matter,new Function::Polynomial(coefficients2),MobilizedBodyIndex(1), MobilizerQIndex(5));
   */


	/*Real value1 = -107.0*Pi/180;
	int argumentSize = 1;
	Constraint::PrescribedMotion constraint2(matter,new Function::Constant(value1, argumentSize),MobilizedBodyIndex(2), MobilizerQIndex(0));
	
	Real value2 = -106.0*Pi/180;
	Constraint::PrescribedMotion constraint3(matter,new Function::Constant(value2, argumentSize),MobilizedBodyIndex(3), MobilizerQIndex(0));*/

	//Constraint::PrescribedMotion constraint4(matter,new Function::Constant(value1, argumentSize),MobilizedBodyIndex(17), MobilizerQIndex(0));

	//Constraint::PrescribedMotion constraint5(matter,new Function::Constant(value1, argumentSize),MobilizedBodyIndex(15), MobilizerQIndex(0));

	//Constraint::PrescribedMotion constraint6(matter,new Function::Constant(value1, argumentSize),MobilizedBodyIndex(14), MobilizerQIndex(0));

	//Constraint::PrescribedMotion constraint7(matter,new Function::Constant(value1, argumentSize),MobilizedBodyIndex(13), MobilizerQIndex(0));

	//Constraint::PrescribedMotion constraint8(matter,new Function::Constant(value1, argumentSize),MobilizedBodyIndex(12), MobilizerQIndex(0));
	
	//Constraint::PrescribedMotion constraint9(matter,new Function::Constant(value1, argumentSize),MobilizedBodyIndex(11), MobilizerQIndex(0));
	
	//Constraint::PrescribedMotion constraint10(matter,new Function::Constant(value1, argumentSize),MobilizedBodyIndex(10), MobilizerQIndex(0));

    //Constraint::PrescribedMotion constraint11(matter,new Function::Constant(value1, argumentSize),MobilizedBodyIndex(9), MobilizerQIndex(0));

	//Constraint::PrescribedMotion constraint12(matter,new Function::Constant(value1, argumentSize),MobilizedBodyIndex(8), MobilizerQIndex(0));

	

    //   Real value1 = -74.6426*Pi/180;;
	//int argumentSize = 1;
	//Constraint::PrescribedMotion constraint2(matter,new Function::Constant(value1, argumentSize),MobilizedBodyIndex(16), MobilizerQIndex(0));
	//
	//Real value2 = -73.4752*Pi/180;;
	//Constraint::PrescribedMotion constraint3(matter,new Function::Constant(value2, argumentSize),MobilizedBodyIndex(14), MobilizerQIndex(0));
	//
	//Real value3 = -70.1764*Pi/180;;
	//Constraint::PrescribedMotion constraint4(matter,new Function::Constant(value3, argumentSize),MobilizedBodyIndex(11), MobilizerQIndex(0));
	//
	//
	//Real value = 0;
	//Constraint::PrescribedMotion constraint5(matter,new Function::Constant(value, argumentSize),MobilizedBodyIndex(7), MobilizerQIndex(0));
 //   //Constraint::PrescribedMotion constraint3(matter,new Function::Constant(value, argumentSize),MobilizedBodyIndex(4), MobilizerQIndex(0));
 //   Constraint::PrescribedMotion constraint6(matter,new Function::Constant(value, argumentSize),MobilizedBodyIndex(5), MobilizerQIndex(0));
	//Constraint::PrescribedMotion constraint7(matter,new Function::Constant(value, argumentSize),MobilizedBodyIndex(12), MobilizerQIndex(0));

	
        state = system.realizeTopology(); // The prescribed motion is now applied on the system
        system.realize(state,Stage::Position); // state::Position is necessary to invoke findBodyTransformInAnotherBody 
        // This part is added to constrain the last residue in the system welded to the ground
		int nBodies;
		nBodies = system.getMatterSubsystem().getNBodies();
        MobilizedBody body_last = matter.getMobilizedBody(MobilizedBodyIndex(nBodies-1));

		Constraint::Weld constraint13(body_last,Transform(Vec3(0)) , 
                   matter.Ground(), 
                   (body_last).findBodyTransformInAnotherBody(state,matter.Ground()) ); // I got the idea from the weldtograound part
		state = system.realizeTopology(); // The constraint is now applied on the system
        
		
		/*cout << constraint0.getConstraintIndex();
	    cout << constraint1.getConstraintIndex();
	    cout << constraint2.getConstraintIndex();
	    cout << constraint3.getConstraintIndex();
		cout << constraint4.getConstraintIndex();*/


    
		if (myParameterReader.verbose) {
        const State & state5 = state; 
        string myFileName5 = "check.5.pdb";
        ofstream  myOfstream5(  myFileName5.c_str());
	for (int i = 0; i < (myParameterReader.sequences).size(); i++) {
            myMolecule[i].writePdb(state5,( myOfstream5),Transform(Vec3(0)));
            if (myParameterReader.verbose) cout<<"[Repel.h] check 5 writing molecule i="<<i<<endl;   
        }
        time ( &rawtime );
        timeinfo = localtime ( &rawtime );
        if (myParameterReader.verbose) cout<<"[Repel.h] check 3: time:"<< asctime (timeinfo)   <<endl;
        }

	Vector qVector;
        qVector=state.getQ();
	Vector oldQVector = qVector;
        if (myParameterReader.verbose) cout<<qVector<<endl<<"The above is the qVector fitted from PDB file"<<endl;
        if (myParameterReader.readInQVector) 
        { 
	    if (myParameterReader.verbose) cout<<"[Repel.h] myParameterReader.inQVectorFileName ="<<myParameterReader.inQVectorFileName<<endl;	
            ifstream lastQVectorFile((myParameterReader.inQVectorFileName).c_str());
            char  * w;
	    string line;
	    int i = 0;
	    while (lastQVectorFile.good()) 
            {
	    	getline(lastQVectorFile,line);//lastQVectorFile.getline(w,100);
		if (myParameterReader.verbose) cout<<"[Repel.h] just read :"<<line <<endl;
		qVector[i] = atof(line.c_str());
	        i++;
	    }
	    lastQVectorFile.close();
            if (myParameterReader.verbose) cout<<qVector<<endl<<"The above is the qVector read from Q vector file"<<endl;
	assert((oldQVector).size() == qVector.size());
        state.setQ(qVector);


        if (myParameterReader.verbose) cout<<state.getQ()<<endl<<"The above is the qVector reset from qVector file"<<endl;
        }
/*
*/
	//myMolecule[0].writePdb(state,ostream('after-realizePosition.pdb'),Vec3(0));

 /*        
        state = system.realizeTopology();
        system.realize(state,Stage::Position);
*/
        if (myParameterReader.weldToGround) { 
            //const SimTK::Transform& mytransform2 = matter.getMobilizedBody(myMolecule[0].getAtomMobilizedBodyIndex(Compound::AtomIndex(0))).findBodyTransformInAnotherBody(state,matter.getMobilizedBody(myMolecule[0].getAtomMobilizedBodyIndex(Compound::AtomIndex(myMolecule[0].getNAtoms()-1))));
            vector<Constraint> myWeld4(myParameterReader.numChains);//[myParameterReader.numRigidSegments];
            for (int i = 0; i<myParameterReader.numChains; i++) {
            if (myParameterReader.verbose) cout<<"[Repel.h] myMolecule[i].getAtomMobilizedBodyIndex(myMolecule[i].getAtomIndex(0/C3*))) ="<<myMolecule[i].getAtomMobilizedBodyIndex(myMolecule[i].getAtomIndex("0/C3*"))<<endl;
            const SimTK::Transform& mytransform3 = (matter.getMobilizedBody(myMolecule[i].getAtomMobilizedBodyIndex(myMolecule[i].getAtomIndex("0/C3*")))).findBodyTransformInAnotherBody(state,matter.Ground());            
                myWeld4[i] =Constraint::Weld(
                    matter.updMobilizedBody(myMolecule[i].getAtomMobilizedBodyIndex(Compound::AtomIndex(myMolecule[i].getAtomIndex("0/C3*") ))),Transform(Vec3(0)) , 
                    matter.Ground(), 
                    (matter.getMobilizedBody(myMolecule[i].getAtomMobilizedBodyIndex(myMolecule[i].getAtomIndex("0/C3*")))).findBodyTransformInAnotherBody(state,matter.Ground())  
                //mytransform3
                );  
                state = system.realizeTopology();
                system.realize(state,Stage::Position);
            }
        }

        // now we will weld together all paired RigidSegments
                                    //const SimTK::Transform mytransform3 = (matter.Ground()).findBodyTransformInAnotherBody(state,(matter.Ground())); // this was ok

        Constraint myWeld1[maxChiWelds];//[myParameterReader.numRigidSegments];
        //                            const SimTK::Transform mytransform3 = (matter.Ground()).findBodyTransformInAnotherBody(state,(matter.Ground())); // this was ok
        if (0) if (myParameterReader.rigidifyFormedHelices) for (int i = 0; i < (myParameterReader.sequences).size(); i++) {
            RNA& rna = static_cast<RNA&>(myMolecule[i]);
            RigidifyFormedHelices(myParameterReader,rna,matter,system,state,
	        (&(myWeld1[0])));
            myMolecule[i]=rna;
        }
        time ( &rawtime );
        timeinfo = localtime ( &rawtime );
        if (myParameterReader.verbose) cout<<"[Repel.h] check 4: time:"<< asctime (timeinfo)   <<endl;
        state = system.realizeTopology();
        system.realize(state,Stage::Position);


        time ( &rawtime );
        timeinfo = localtime ( &rawtime );
        if (myParameterReader.verbose) cout<<"[Repel.h] check 5: time:"<< asctime (timeinfo)   <<endl;
	
         //                           const SimTK::Transform mytransform3 = (matter.Ground()).findBodyTransformInAnotherBody(state,(matter.Ground())); // this was not ok


        //vector<Biopolymer> myMoleculeVector(myMolecule, myMolecule  + sizeof(myMolecule ) / sizeof(Biopolymer) );
	//for (int i = 0; i<maxChains; i++) myMoleculeVector[i] = myMolecule[i];
        vector<Constraint> myWeld2(maxChiWelds);//[myParameterReader.numRigidSegments];
        //RigidifyFormedHelices(myParameterReader,myMoleculeVector,matter,state,myWeld2);
	int r = 0;
        if (1) for (int q=0;q<myParameterReader.numBasePairs;q++) 
        if ((myParameterReader.firstBPEdge[q]).compare("Weld" ) == 0) 
           {    
            int i = -11111;
            int j = -11111;
            for (int k = 0; k < maxChains; k++) {
                if (strcmp((myParameterReader.chainId[k]).c_str(),(myParameterReader.firstBPChain[q]).c_str()) == 0) i = k; 
                if (strcmp((myParameterReader.chainId[k]).c_str(),(myParameterReader.secondBPChain[q]).c_str()) == 0) j = k; 
            }    
            assert(i>=0);
            assert(j>=0);
                                {   
                                    //weld a residue from r with a residue from q
                                    stringstream ss3;
                                    ss3<<myParameterReader.firstBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()) <<"/"<<"C3*";
                                    stringstream ss4;
                                    ss4<<myParameterReader.secondBPResidue[q]-(myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[j]).c_str()))<<"/"<<"C3*";    
				    if (myParameterReader.verbose) cout<<"[Repel.h] ss3,ss4,i,j ="<<ss3.str()<<","<<ss4.str()<<","<<i<<","<<j<<endl;
				    if (myParameterReader.verbose) cout<<(matter.updMobilizedBody((myMolecule[i]).getAtomMobilizedBodyIndex(myMolecule[i].getAtomIndex(ss3.str())))).getMobilizedBodyIndex()<<endl;
				    if (myParameterReader.verbose) cout<<(matter.updMobilizedBody((myMolecule[j]).getAtomMobilizedBodyIndex(myMolecule[j].getAtomIndex(ss4.str())))).getMobilizedBodyIndex()<<endl;
                                    //const SimTK::Transform mytransform2 = (matter.updMobilizedBody((myMolecule[i]).getAtomMobilizedBodyIndex((myMolecule[i]).getAtomIndex(ss3.str())))).findBodyTransformInAnotherBody(state,(matter.updMobilizedBody((myMolecule[j]).getAtomMobilizedBodyIndex((myMolecule[j]).getAtomIndex(ss4.str())))));
                                    const SimTK::Transform mytransform1 = (matter.updMobilizedBody((myMolecule[i]).getAtomMobilizedBodyIndex((myMolecule[i]).getAtomIndex(ss3.str())))).findBodyTransformInAnotherBody(state,(matter.updMobilizedBody((myMolecule[j]).getAtomMobilizedBodyIndex((myMolecule[j]).getAtomIndex(ss4.str())))));
    //Finally, we must Weld the two RigidSegments together, since the forces will no longer work.
				    assert(r<maxChiWelds);
				    if (myParameterReader.verbose) cout<<"[Repel.h] r="<<r<<endl;
                                    myWeld2[r] = Constraint::Weld  (
                                        matter.updMobilizedBody((myMolecule[i]).getAtomMobilizedBodyIndex((myMolecule[i]).getAtomIndex(ss3.str()))),
                                        Transform(Vec3(0)),
                                        (matter.updMobilizedBody((myMolecule[j]).getAtomMobilizedBodyIndex((myMolecule[j]).getAtomIndex(ss4.str())))),
                                        mytransform1
                                    );  
				    r++;
                                    state = system.realizeTopology();
                                    system.realize(state,Stage::Position);

                                } 

           }

        AllTwoTransformLinearSprings * myAllTwoTransformLinearSpringsPointer = 
        new AllTwoTransformLinearSprings( matter,  myParameterReader,  myLeontisWesthofClass, &(myMolecule[0]));
        //Force::Custom(forces, myAllTwoTransformLinearSpringsPointer);
            //discontinued setChiBondAnti, since it's the only thing that uses the now-discontinued BasePairingBonds.
            /*
            if (0) if (myParameterReader.setChiBondAnti) for (int i = 0; i < (myParameterReader.sequences).size(); i++) { 
            for (int q=0;q<myParameterReader.numRigidSegments[(myParameterReader.chainId[i]).c_str()]; q++){ 
                if (myParameterReader.BPEdge[(myParameterReader.rigidSegmentStartNumber[q][i]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()))][i].compare("WatsonCrick") ==0 ){
		    for (int myResidueNumber =myParameterReader.rigidSegmentStartNumber[q][i]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()) ;myResidueNumber <=  myParameterReader.rigidSegmentEndNumber[q][i]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()); myResidueNumber++) {
			//cout<<"[Repel.h] myParameterReader.rigidSegmentStartNumber[q][i],  myParameterReader.rigidSegmentEndNumber[q][i] ="<< myParameterReader.rigidSegmentStartNumber[q]<<","  << myParameterReader.rigidSegmentEndNumber[q]<<endl;
                        BasePairingBonds(myMolecule,myResidueNumber,myResidueNumber,myParameterReader.chainId[i],myParameterReader.chainId[i],"ChiBondAnti","ChiBondAnti","Cis",myLeontisWesthofClass,forces,matter,myHuntCrossleyContact,contacts, contactSet,hc ,state,myParameterReader.dutyCycle , myParameterReader);
		    }   
                }
            }
            }*/
	//this is a force that gets all the base pairing info it needs from myParameterReader and myLeontisWesthofClass.  This way myParameterReader can be updated and the forces will immediately be updated.


	// this polymorphism has its own loop to count over base pairs:
        //BasePairingBonds(myMolecule,myLeontisWesthofClass,forces,matter,myHuntCrossleyContact,contacts, contactSet,hc ,state,myParameterReader.dutyCycle,myParameterReader);
        /*Force::TwoPointLinearSpring myTwoPointLinearSpring;

        if (myParameterReader.addTestSpring) myTwoPointLinearSpring= Force::TwoPointLinearSpring(
            forces,
            matter.updMobilizedBody(myMagnesiumIonVec[0].getAtomMobilizedBodyIndex(myMagnesiumIonVec[0].getAtomIndex("MG+"))),
            Vec3(0),
            matter.updMobilizedBody(myMagnesiumIonVec[1].getAtomMobilizedBodyIndex(myMagnesiumIonVec[1].getAtomIndex("MG+"))),
            Vec3(0) ,
            1.0,
	    0.0
            
            );
*/
        if (0) if (myParameterReader.addTestSpring)  
        Force::TwoPointLinearSpring myTwoPointLinearSpring(
            forces,
            matter.updMobilizedBody(myMagnesiumIonVec[0].getAtomMobilizedBodyIndex(myMagnesiumIonVec[0].getAtomIndex("MG+"))),
            Vec3(0),
            matter.updMobilizedBody(myMagnesiumIonVec[1].getAtomMobilizedBodyIndex(myMagnesiumIonVec[1].getAtomIndex("MG+"))),
            Vec3(0) ,
            1.0,
	    0.0
            
            );
        Force::Custom(forces, myAllTwoTransformLinearSpringsPointer);
//	if (0) 
        for (int q=0;q<myParameterReader.numBasePairs;q++) 
        if  //((myParameterReader.basePairIsTwoTransformForce[q]).compare("contact") == 0)   
		(((myParameterReader.firstBPEdge[q]).compare("HardSphere") == 0) ||
		 ((myParameterReader.firstBPEdge[q]).compare("HardSphereSmall") == 0) ||
		 ((myParameterReader.firstBPEdge[q]).compare("AddAllHeavyAtomSterics") == 0) ||
		 ((myParameterReader.firstBPEdge[q]).compare("ExcludedVolume") == 0) )
           {
		            int i = -11111;
            int j = -11111;
            for (int k = 0; k < maxChains; k++) {
                if (strcmp((myParameterReader.chainId[k]).c_str(),(myParameterReader.firstBPChain[q]).c_str()) == 0) i = k;
                if (strcmp((myParameterReader.chainId[k]).c_str(),(myParameterReader.secondBPChain[q]).c_str()) == 0) j = k;
            }   
            assert(i>=0);
            assert(j>=0);
	    //assert (i == j);
                SimTK_ERRCHK_ALWAYS(
                    ((i == j)) ,
                    "[Repel.h]",
                    "For HardSphere and HardSphereSmall, you must have the start and end points on the same chain."); 

            assert ((myParameterReader.firstBPEdge[q]).compare(myParameterReader.firstBPEdge[q]) == 0);

		if (myParameterReader.verbose) cout<<"about to do BasePairingBonds :"<<myParameterReader.firstBPResidue[q]<<","<<myParameterReader.secondBPResidue[q]<<endl;
		if (myParameterReader.verbose) cout<<"[Repel.h] myParameterReader.firstBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()),myParameterReader.secondBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[j]).c_str()) ="<<myParameterReader.firstBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str())<<","<<myParameterReader.secondBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[j]).c_str())<<endl;
		if (myParameterReader.verbose) cout<<"[Repel.h] myParameterReader.firstBPResidue[q],myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()),myParameterReader.secondBPResidue[q],myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[j]).c_str()) ="<<myParameterReader.firstBPResidue[q]<<","<<myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str())<<","<<myParameterReader.secondBPResidue[q]<<","<<myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[j]).c_str())<<endl;
                if ((myParameterReader.firstBPEdge[q]).compare("ExcludedVolume") ==0 )
                    AddExcludedVolume     (myMolecule[i],contacts,contactSetLargeSpheres ,hcLargeSpheres, matter,myParameterReader    ,myParameterReader.firstBPEdge[q], myParameterReader.firstBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()),myParameterReader.secondBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()));
                else if ((myParameterReader.firstBPEdge[q]).compare("AllHeavyAtomSterics") ==0 )
                    AddAllHeavyAtomSterics     (myMolecule[i],contacts,contactSetLargeSpheres ,hcLargeSpheres, matter,myParameterReader    ,myParameterReader.firstBPEdge[q], myParameterReader.firstBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()),myParameterReader.secondBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()));
                     
                else AddHuntCrossleySpheres(myMolecule[i],contacts,contactSetLargeSpheres ,hcLargeSpheres, matter,myLeontisWesthofClass,myParameterReader.firstBPEdge[q], myParameterReader.firstBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()),myParameterReader.secondBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()));
                    //AddHuntCrossleySpheres(myMolecule[0],contacts,contactSetLargeSpheres,hcLargeSpheres , matter,myLeontisWesthofClass,"HardSphere");
               // BasePairingBonds(myMolecule, myParameterReader.firstBPResidue[q], myParameterReader.secondBPResidue[q], myParameterReader.chainId[i],myParameterReader.chainId[j] , myParameterReader.firstBPEdge[q], myParameterReader.secondBPEdge[q], myParameterReader.orientationBP[q], myLeontisWesthofClass, forces, matter, myHuntCrossleyContact, contacts, contactSet, hc, state, myParameterReader.dutyCycle, myParameterReader);

           }
           else if ((myParameterReader.firstBPEdge[q]).compare("BackboneOxygen") == 0){
	       int i = -11111;
               for (int k = 0; k < maxChains; k++) if (strcmp((myParameterReader.chainId[k]).c_str(),(myParameterReader.firstBPChain[q]).c_str()) == 0) i = k;
               myAddBackboneOxygenForces.addForces( system, matter , backboneOxygenForces, myParameterReader, myMolecule[i],myParameterReader.firstBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()),myParameterReader.secondBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()) );
           }
        //for (int i = 0; i<maxChains; i++) myMolecule[i] = myMoleculeVector[i];
        int i =0;
        for (int q=0;q<myParameterReader.numBasePairs;q++) 
        if ((myParameterReader.firstBPEdge[q]).compare("HardSphere") == 0)  {
            for (int k = 0; k < maxChains; k++) 
                if (strcmp((myParameterReader.chainId[k]).c_str(),(myParameterReader.firstBPChain[q]).c_str()) == 0) i = k;
                AddHuntCrossleySpheres(myMolecule[i],contacts,contactSetLargeSpheres ,hcLargeSpheres, matter,myLeontisWesthofClass,myParameterReader.firstBPEdge[q], myParameterReader.firstBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()),myParameterReader.secondBPResidue[q]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str()));
        }


 
       if (myParameterReader.addExcludedVolume     ) {
		for (int m = 0; m < (myParameterReader.sequences).size(); m++) 
                    {AddExcludedVolume     (myMolecule[m],contacts,contactSetLargeSpheres,hcLargeSpheres , matter,myParameterReader    ,"ExcludedVolume",0,myMolecule[m].getNumResidues() -1 );
                    cout <<"";
                    }
       }

       if (myParameterReader.addHuntCrossleySpheres) {
		for (int m = 0; m < (myParameterReader.sequences).size(); m++) 
                    {AddHuntCrossleySpheres(myMolecule[m],contacts,contactSetLargeSpheres,hcLargeSpheres , matter,myLeontisWesthofClass,"HardSphere");
                    cout <<"";
                    }
            
            for (int i = 0 ; i< myParameterReader.numMagnesiumIons; i++) {
                contacts.addBody(contactSetLargeSpheres,
                                        (matter.updMobilizedBody(myMagnesiumIonVec[i].getAtomMobilizedBodyIndex(Compound::AtomIndex(0)))), 
                                        ContactGeometry::Sphere(.1), 
                                        Vec3(0) //(myMagnesiumIonVec[i].getAtomLocationInMobilizedBodyFrame(myMolecule.getAtomIndex(ss3.str())))
                                        );  
                hcLargeSpheres.setBodyParameters(contacts.getNumBodies(contactSetLargeSpheres)-1,10000,0, 0., 0., 0.);
            }
            ContactGeometry::HalfSpace maxXHalfSpace;
        /*    HalfSpace minXHalfSpace;
            HalfSpace maxYHalfSpace;
            HalfSpace minYHalfSpace;
            HalfSpace maxZHalfSpace;
            HalfSpace minZHalfSpace;*/
            MobilizedBody myMobilizedBody[6];
            contacts.addBody(contactSetLargeSpheres,
                matter.Ground(),   
                ContactGeometry::HalfSpace(),
                Transform(Rotation(0,Vec3(0,0,1)),Vec3(myParameterReader.halfSpaceMaxX,0,0))
                ); 
            hcLargeSpheres.setBodyParameters(contacts.getNumBodies(contactSetLargeSpheres)-1,10000,0, 0., 0., 0.);
            contacts.addBody(contactSetLargeSpheres,
                matter.Ground(),   
                ContactGeometry::HalfSpace(),
                Transform(Rotation(180/Rad2Deg,Vec3(0,0,1)),Vec3(-myParameterReader.halfSpaceMaxX,0,0))
                );  
            hcLargeSpheres.setBodyParameters(contacts.getNumBodies(contactSetLargeSpheres)-1,10000,0, 0., 0., 0.);
            contacts.addBody(contactSetLargeSpheres,
                matter.Ground(),   
                ContactGeometry::HalfSpace(),
                Transform(Rotation(90/Rad2Deg,Vec3(0,0,1)),Vec3(0,myParameterReader.halfSpaceMaxX,0))
                );  
            hcLargeSpheres.setBodyParameters(contacts.getNumBodies(contactSetLargeSpheres)-1,10000,0, 0., 0., 0.);
            contacts.addBody(contactSetLargeSpheres,
                matter.Ground(),   
                ContactGeometry::HalfSpace(),
                Transform(Rotation(-90/Rad2Deg,Vec3(0,0,1)),Vec3(0,-myParameterReader.halfSpaceMaxX,0))
                );  
            hcLargeSpheres.setBodyParameters(contacts.getNumBodies(contactSetLargeSpheres)-1,10000,0, 0., 0., 0.);
            contacts.addBody(contactSetLargeSpheres,
                matter.Ground(),   
                ContactGeometry::HalfSpace(),
                Transform(Rotation(-90/Rad2Deg,Vec3(0,1,0)),Vec3(0,0,myParameterReader.halfSpaceMaxX))
                );  
            hcLargeSpheres.setBodyParameters(contacts.getNumBodies(contactSetLargeSpheres)-1,10000,0, 0., 0., 0.);
            contacts.addBody(contactSetLargeSpheres,
                matter.Ground(),   
                ContactGeometry::HalfSpace(),
                Transform(Rotation( 90/Rad2Deg,Vec3(0,1,0)),Vec3(0,0,-myParameterReader.halfSpaceMaxX))
                ); 
            hcLargeSpheres.setBodyParameters(contacts.getNumBodies(contactSetLargeSpheres)-1,10000,0, 0., 0., 0.);
/*
 */
            
        }
        if (myParameterReader.addSmallHuntCrossleySpheres) 
		for (int i = 0; i < (myParameterReader.sequences).size() ; i++ ) 
		    AddHuntCrossleySpheres(myMolecule[i],contacts,contactSet,hc, matter,myLeontisWesthofClass,"HardSphereSmall");
//AddSmallHuntCrossleySpheres(myMolecule[i],contacts,contactSet,hc, matter,myLeontisWesthofClass,"HardSphereSmall");


        if (myParameterReader.setTemperature){

        // last param is relax time in ps.  if this number is small, integrator will be slow.  this is the time it takes to enforce temperature.
        // new nose-hoover.  uncomment this when Simbody has been updated. then get rid of velocityrescaling thermostat
        //
        if (myParameterReader.verbose) cout<<"[Repel.h] myParameterReader.temperature *100="<<myParameterReader.temperature*100<<endl;
        if ((myParameterReader.thermostatType).compare("NoseHoover") == 0 ) 
            Force::Thermostat  (forces, matter,SimTK_BOLTZMANN_CONSTANT_MD,myParameterReader.temperature,1.0);
        else if ((myParameterReader.thermostatType).compare("VelocityRescaling") == 0 ) { 
            VelocityRescalingThermostat * myVelocityRescalingThermostat = new VelocityRescalingThermostat(system,myParameterReader.temperature, myParameterReader.velocityRescalingInterval);
            system.updDefaultSubsystem().addEventHandler(myVelocityRescalingThermostat);
        }
            }
            //VelocityRescalingThermostat * myVelocityRescalingThermostat = new VelocityRescalingThermostat(system,    temperature, myReportingInterval);
        else
            Force::GlobalDamper(forces, matter,.3);
        ofstream output(myParameterReader.outTrajectoryFileName.c_str());//"mymovie.2.pdb");
        PeriodicPdbAndEnergyWriter * myPeriodicPdbWriter = new PeriodicPdbAndEnergyWriter(system,output,  myParameterReader.reportingInterval);
        PeriodicParameterReaderUpdater * myPeriodicParameterReaderUpdater = new PeriodicParameterReaderUpdater(system,*myAllTwoTransformLinearSpringsPointer ,  myParameterReader.reportingInterval*10);
	char * frameFileName = "frame.pdb";
        /*
	char * frameFileName = new char[myParameterReader.lastFrameFileName.size()+1] ;
        (myParameterReader.lastFrameFileName).copy(frameFileName,(myParameterReader.lastFrameFileName).size(),0) ;
        frameFileName[myParameterReader.lastFrameFileName.size()] = 0;
        */
        PeriodicPdbAndEnergySingleFrameWriter * myPeriodicPdbSingleFrameWriter = new PeriodicPdbAndEnergySingleFrameWriter(system,frameFileName  ,  myParameterReader.reportingInterval);
        
        //PeriodicPdbWriter * myPeriodicPdbWriter = new PeriodicPdbWriter(system,output,  myParameterReader.reportingInterval);
        PeriodicScrubber  * myPeriodicScrubber  = new PeriodicScrubber (system,forces,  0.4*10,myParameterReader.dutyCycle);
        //VTKEventReporter * myVTKEventReporter = new VTKEventReporter(system,myParameterReader.reportingInterval);
        //VTKEventReporter * myVTKEventReporter = new VTKEventReporter(system,myParameterReader.reportingInterval);

        //Event handler to relase the joint at the specific time
         system.updDefaultSubsystem().addEventHandler(new release_joint(matter));
		// state = system.realizeTopology();

		//Event handler to lock the joint at the specific time
         //system.updDefaultSubsystem().addEventHandler(new lock_joint(matter));

       for (int q=0; q<myParameterReader.numBasePairs; q++) if (myParameterReader.verbose) cout<<"Base pair number, firstBPResidue, firstBPEdge,secondBPResidue,secondBPEdge, priority ="<<q <<","<< + myParameterReader.firstBPResidue[q]<<","<< myParameterReader.firstBPEdge[q] <<","<<myParameterReader.secondBPResidue[q]<<","<<myParameterReader.secondBPEdge[q]<<","<<   myParameterReader.basePairPriority[q]<<endl; 	
        system.updDefaultSubsystem().addEventReporter(myPeriodicPdbWriter);
        system.updDefaultSubsystem().addEventReporter(myPeriodicPdbSingleFrameWriter);
        if (myParameterReader.periodicallyUpdateParameters) system.updDefaultSubsystem().addEventHandler (myPeriodicParameterReaderUpdater);
        if (myParameterReader.setForceAndStericScrubber) system.updDefaultSubsystem().addEventHandler (myPeriodicScrubber );
        //system.updDefaultSubsystem().addEventReporter(new VTKEventReporter(system,myParameterReader.reportingInterval));//myVTKEventReporter);
            
       // Send coordinates to VMD 
        if (myParameterReader.vmdOutput == 1) 
            system.updDefaultSubsystem().addEventReporter(new PeriodicVmdReporter(
                system, 0.015, 3000, true));

 
        //CentrifugalForce myCentrifugalForce(matter);
        //StackingTorque  myStackingTorque(matter,myMolecule[0]);
        /*
        Integrator study;
        if (useRunge) {RungeKuttaMersonIntegrator study(system);//,myParameterReader.stepSize);
        VerletIntegrator study(system);//,myParameterReader.stepSize);
        //VerletIntegrator study;
        if (myParameterReader.useFixedStepSize)
        //else 
          //  study = VerletIntegrator(system);//,myParameterReader.stepSize);
            //study = RungeKuttaMersonIntegrator(system);*/
        if (myParameterReader.monteCarloRun) {
            cout<<"[Repel.h] About to start Monte Carlo RandomizeRNACoordinates run "<<endl;
            RandomizeRNACoordinates  ( system,  state, matter ,myParameterReader  ,static_cast<RNA&>(myMolecule[0]));
        }
        Integrator *study;
        if ((myParameterReader.integratorType).compare("RungeKuttaMerson")) 
            study = new RungeKuttaMersonIntegrator(system);
        else 
            if ((myParameterReader.integratorType).compare("Verlet"          ))
                study = new VerletIntegrator(system);
            else
                {cout<<"[Repel.h] You must specify an integratorType of one of the supported types!!!"<<endl;assert(0);}
        if (myParameterReader.integratorUseFixedStepSize)
            study->setFixedStepSize(myParameterReader.integratorStepSize);
        study->setAccuracy(myParameterReader.integratorAccuracy);
        study->setConstraintTolerance(myParameterReader.constraintTolerance);//myParameterReader.integratorAccuracy);
        state = system.realizeTopology();
        /*
        for (int i = 0; i < forces.getNForces(); i++) { 
            forces.setForceIsDisabled(state,SimTK::ForceIndex(i),1);
	}*/
	if (myParameterReader.minimize) {
            cout<<"about to start minimizing"<<endl; 

            if (myParameterReader.verbose) cout<<state.getQ()<<endl<<"The above is the qVector before resetting                              "<<endl;
            for (int k = 0; k<6; k++) state.updQ()[k]=1+k;//Vector(0.,1.,2.,3.,4.,5.));
            if (myParameterReader.verbose) cout<<state.getQ()<<endl<<"The above is the qVector before energy minimization                    "<<endl;
            LocalEnergyMinimizer::minimizeEnergy(system,state, .1); 
            if (myParameterReader.verbose) cout<<state.getQ()<<endl<<"The above is the qVector after energy minimization e"<<endl;
            //assert(0);
	    if (myParameterReader.verbose) cout<<"done minimizing "<<endl;
            Vector_<SpatialVec>  myBodyForceVector(3); 
            Vector_<Vec3      >  myBodyParticleForceVector(3); 
            Vector           	 myMobilityForceVector(3); 
            //(myTwoPointLinearSpring.getImpl()).calcForce(state,myBodyForceVector,myBodyParticleForceVector,myMobilityForceVector  );
            if (myParameterReader.verbose) cout<<"myBodyForceVector         = "<<myBodyForceVector        <<endl;
            if (myParameterReader.verbose) cout<<"myMobilityForceVector = "<<myMobilityForceVector<<endl;
        }
	if (myParameterReader.verbose) cout<<state.getQ()<<"The above is qVector"<<endl;
        if (0) if (myParameterReader.readInQVector) {
	    state.setQ(qVector);
            if (myParameterReader.verbose) cout<<state.getQ()<<endl<<"The above is the qVector which shoudl have been reset from qVector file"<<endl;
 	}    


        time ( &rawtime );
        timeinfo = localtime ( &rawtime );
        if (myParameterReader.verbose) cout<<"[Repel.h] check 6: time:"<< asctime (timeinfo)   <<endl;
       // state = system.realizeTopology();
        //system.realize(state,Stage::Position);
        TimeStepper ts(system,*study);
        ts.initialize(state);
        cout<<"[Repel.h:ConstrainedDynamics] Starting dynamics now."<<endl;
        ts.stepTo(myParameterReader.maxReportingIntervals*myParameterReader.reportingInterval);
        //delete study; //peter
/*
        state=   ts.getState();
        string myFileName13 = "check.13.pdb";
        ofstream  myOfstream13(  myFileName13.c_str());
	for (int i = 0; i < (myParameterReader.sequences).size(); i++) {
            myMolecule[i].writePdb(state,( myOfstream13),Transform(Vec3(0)));
            cout<<"[Repel.h] check 13 writing molecule i="<<i<<endl;   
        }
*/
  /*      cout<<"Crrently using constraint tolerance ="<<study->getConstraintToleranceInUse()<<endl;
        cout<<"getNStepsAttempted  ="<<study->getNStepsAttempted()<<endl;
        cout<<"getNStepsTaken  ="<<study->getNStepsTaken()<<endl;
        cout<<"Currently using accuracy ="<<study->getAccuracyInUse()<<endl;
	cout<<"Last step size ="<<study->getPreviousStepSizeTaken() <<endl;
    */    
        //if (myParameterReader.readInQVector) 
        {
            state=   ts.getState();
            qVector=state.getQ();//study->getState()).getQ();
	    ofstream lastFrameQVectorFile((myParameterReader.outQVectorFileName).c_str()); 
	    for (int q = 0; q<qVector.size(); q++) lastFrameQVectorFile<<qVector[q]<<endl;
            lastFrameQVectorFile.close();
        }
        /*
	*/
        ofstream lastFrameFile(myParameterReader.lastFrameFileName.c_str()); 
                for (SimTK::Compound::Index c(0); c < system.getNumCompounds(); ++c)
                    (system.getCompound(c)).writePdb(state, lastFrameFile,Transform(Vec3(0)));//, nextAtomSerialNumber);


 //       for (int i = 0; i < (myParameterReader.sequences).size(); i++) myMolecule[i].writePdb(state,lastFrameFile,Transform(Vec3(0)));
       
        lastFrameFile.close();
        //stringstream ss9;
	//ss9.clear();
        //ss9<<"awk \'{if (substr($3,0,1) != \"H\") {print $0}}\' "<<myParameterReader.lastFrameFileName<<" > temp.zzz.pdb ; mv temp.zzz.pdb  "<<myParameterReader.lastFrameFileName<<endl;
        //syscall(ss9); 
        
        ifstream myLastFrameFile(myParameterReader.lastFrameFileName.c_str()); 
        ofstream lastFrameFileTemp("temp.zzz.pdb"           ); 
        char * v;
        v = new char[1000] ;
        stringstream u;
        string mystring[10];
        while (myLastFrameFile.good()) {
            myLastFrameFile.getline(v,1000);
            //cout <<"[ParameterReader.h] v="<<v<<endl;  
            u.clear();
            u.str("");
            u<<string(v);
            u>>mystring[0]>>mystring[1]>>mystring[2]>>mystring[3]>>mystring[4]>>mystring[5]>>mystring[6]>>mystring[7]>>mystring[8]>>mystring[9]; 
            //if ((((mystring[2]).substr(0,1)).compare("H") != 0) && (((mystring[2]).substr(0,1)).compare("N9") != 0) && (((mystring[2]).substr(0,1)).compare("N1") != 0) ) {lastFrameFileTemp<<v<<endl;}
            //cout <<"[ParameterReader.h] myParameterReader.resetBases["<<mystring[4] <<"]="<<myParameterReader.resetBases[mystring[4]]<<endl;  
            
            if (myParameterReader.resetBases[mystring[4]]) 
            {
		if ((((mystring[2]).substr(0,1)).compare("H") != 0) && (((mystring[2])).compare("N9") != 0) && (((mystring[2])).compare("N1") != 0) ) {lastFrameFileTemp<<v<<endl;}
		
            } else  {
                //if ((((mystring[2]).substr(0,1)).compare("H") != 0) ) 
                    {lastFrameFileTemp<<v<<endl;}
            }

        }
        lastFrameFileTemp.close();
        std::ifstream ifs("temp.zzz.pdb", std::ios::binary);
        //std::ofstream ofs(myParameterReader.lastFrameFileName.c_str(), std::ios::binary);
        //ofs << ifs.rdbuf(); 
        //ofs.close();
        output.close();
	return myMolecule[0];
    };	
};   


int PushBasesApart (char * inFileName="assembled.singlechain.pdb") {
    int firstres=17 ;
    float springConstant = 100;
    float deadLength     = .3 ;
    CompoundSystem system;
    SimbodyMatterSubsystem  matter(system);
    DuMMForceFieldSubsystem dumm(system);
    //char* inFileName="1UUI.simple.pdb";
    //char* inFileName="temp.pdb";//"assembled.pdb";//myduplex-joined.2.pdb";
    char* myBondMobilityType;
    float myRigid = 1;
    float myBondStretchGlobalScaleFactor=0.005;
    float myCoulombGlobalScaleFactor = 0.0;
    float myGbsaGlobalScaleFactor = 0;
    float myVdwGlobalScaleFactor = 0.0;
    float myReportingInterval    =0.1 ;
    float maxTime                =myReportingInterval*400  ;//ps
    float thermostatTimeConstant = 0.0;
    int myUseMultithreadedComputation =1; 
    float myTemperature =  290;
    float waterRadius   = 1.00;
    int numWater         =0   ;
    int numDivalents     = 5  ;
    int makeWaterDroplet=1;	

    dumm.loadAmber99Parameters();
    RNA molA      ("GGCAGAUCUGAGCCUGGGAGCUCUCUGCC",0);

          molA.assignBiotypes();
    molA.setPdbChainId('B');
    	
    molA.renumberPdbResidues(17); 





    std::ifstream inFileStream(inFileName,ifstream::in);
    PdbStructure pdbStructure(inFileStream);
    Compound::AtomTargetLocations atomTargetsA = molA.createAtomTargets(pdbStructure);
    //std::cout<<"atomtargest.szie "<<atomTargets.size()<<"versus atoms in myMolecule = "<<myMolecule.getNAtoms()<<std::endl;    
    // Four steps to a perfect match
    cout<<"about to matchDefaultConfiguration for molA"<<endl;
    molA.matchDefaultConfiguration(atomTargetsA);


    molA.setCompoundBondMobility(BondMobility::Free);
    molA.setRNABondMobility(BondMobility::Rigid,17-17,22-17);
    molA.setRNABondMobility(BondMobility::Rigid,23-17,25-17);
    molA.setRNABondMobility(BondMobility::Rigid,26-17,29-17);
    molA.setRNABondMobility(BondMobility::Rigid,30-17,35-17);
    molA.setRNABondMobility(BondMobility::Rigid,36-17,39-17);
    molA.setRNABondMobility(BondMobility::Rigid,40-17,45-17);

    cout<<"about to adopt molA"<<endl;
    system.adoptCompound(molA);

	
    GeneralForceSubsystem forces(system); //this is needed only for the VanderWallSphere in WaterDroplet
    HuntCrossleyContact myHuntCrossleyContact(system);
    Force::Custom(forces, new CentrifugalForce(matter));
    system.modelCompounds();        
    State state = system.realizeTopology();
    system.realize(state,Stage::Position);

    

    //DuplexMaker(molA,17-17,22-17,40-17,45-17 ,forces, matter);
    //DuplexMaker(molA,26-17,29-17,36-17,39-17 ,forces, matter);
    //DuplexMaker(molC, molD,forces, matter);
// start adding spheres

    //AddHuntCrossleySpheres(molA,myHuntCrossleyContact, matter);
// done with spheres

 cout<<"check 2.0"<<endl;
     if (! myRigid)
     {
    	//const SimTK::Transform& mytransform2 = matter.getMobilizedBody(molA.getAtomMobilizedBodyIndex(Compound::AtomIndex(0))).findBodyTransformInAnotherBody(state,matter.getMobilizedBody(molA.getAtomMobilizedBodyIndex(Compound::AtomIndex(molA.getNAtoms()-1))));
    	cout<<"check 2.5"<<endl;
    	//Constraint::Weld myConstraint(matter.updMobilizedBody(molA.getAtomMobilizedBodyIndex(Compound::AtomIndex(0))),Transform(Vec3(0)) , matter.updMobilizedBody(molA.getAtomMobilizedBodyIndex(Compound::AtomIndex(molA.getNAtoms()-1))),mytransform2      );
    }
    cout<<"check 3.0"<<endl;

#ifdef CREATE_VTK_WINDOW
    VTKVisualizer display(system, 0.1);
#endif
        VelocityRescalingThermostat * myVelocityRescalingThermostat = new VelocityRescalingThermostat(system,  myTemperature, myReportingInterval);

    ofstream output("mymovie.pdb");
    PeriodicPdbWriter * myPeriodicPdbWriter = new PeriodicPdbWriter(system,output,myReportingInterval);
    cout<<"check 4.0"<<endl;
    system.updDefaultSubsystem().addEventHandler(myVelocityRescalingThermostat);
    system.updDefaultSubsystem().addEventReporter(myPeriodicPdbWriter);
    cout<<"check 5.0"<<endl;


    dumm.setCoulombGlobalScaleFactor(myCoulombGlobalScaleFactor);	
    dumm.setBondStretchGlobalScaleFactor(myBondStretchGlobalScaleFactor) ;
    dumm.setGbsaGlobalScaleFactor(myGbsaGlobalScaleFactor);	
    dumm.setVdwGlobalScaleFactor(myVdwGlobalScaleFactor);	
    dumm.setUseMultithreadedComputation(myUseMultithreadedComputation);
    Force::Custom(forces, new CentrifugalForce(matter));

    RungeKuttaMersonIntegrator study(system);
    cout<<"check 6.0"<<endl;


    //scf put back in later
    system.realizeTopology();

    int p = 0;
    stringstream ss1;
    ss1<<"out."<<p<<".pdb";
    ofstream outputFrame(ss1.str().c_str());

    cout<<"p = "<<p<<endl;
    for (int i = 0; i< system.getNumCompounds(); i++)
		{
		cout<<"printing out Compound: "<<i<<endl;	
		//system.updCompound(Compound::Index(i)).writePdb(study.getState(),outputFrame,Vec3(0));
		}
    cout<<"done printing out minimized structure"<<endl;	
    cout<<"multithreading set to : "<<dumm.getUseMultithreadedComputation()<<endl;
    cout<<"Gbsa set to : "<<myGbsaGlobalScaleFactor      <<endl;
#ifdef CREATE_VTK_WINDOW
                                        
    display.report(study.getState());
#endif

    TimeStepper ts(system,study);
    ts.initialize(system.getDefaultState());
    ts.stepTo(maxTime);
    return(0);
}



int HealBonds      (char * inFileName="assembled.repelled.singlechain.pdb") {
    int firstres=17 ;
    float springConstant = 100;
    float deadLength     = .3 ;
    CompoundSystem system;
    SimbodyMatterSubsystem  matter(system);
    DuMMForceFieldSubsystem dumm(system);
    //char* inFileName="1UUI.simple.pdb";
    //char* inFileName="temp.pdb";//"assembled.pdb";//myduplex-joined.2.pdb";
    char* myBondMobilityType;
    float myRigid = 1;
    float myBondStretchGlobalScaleFactor=0.05 ;
    float myCoulombGlobalScaleFactor = 0.0;
    float myGbsaGlobalScaleFactor = 0;
    float myVdwGlobalScaleFactor = 0.5;
    float myReportingInterval    =0.1 ;
    float maxTime                =myReportingInterval*400  ;//ps
    float thermostatTimeConstant = 0.0;
    int myUseMultithreadedComputation =1; 
    float myTemperature =  290;
    float waterRadius   = 1.00;
    int numWater         =0   ;
    int numDivalents     = 0  ;
    int makeWaterDroplet=0;	

    dumm.loadAmber99Parameters();
    RNA molA      ("GGCAGAUCUGAGCCUGGGAGCUCUCUGCC",0);

          molA.assignBiotypes();
    molA.setPdbChainId('B');
    	
    molA.renumberPdbResidues(17); 





    std::ifstream inFileStream(inFileName,ifstream::in);
    PdbStructure pdbStructure(inFileStream);
    Compound::AtomTargetLocations atomTargetsA = molA.createAtomTargets(pdbStructure);
    //std::cout<<"atomtargest.szie "<<atomTargets.size()<<"versus atoms in myMolecule = "<<myMolecule.getNAtoms()<<std::endl;    
    // Four steps to a perfect match
    cout<<"about to matchDefaultConfiguration for molA"<<endl;
    molA.matchDefaultConfiguration(atomTargetsA);


    molA.setCompoundBondMobility(BondMobility::Free);
    molA.setRNABondMobility(BondMobility::Rigid,17-17,22-17);
    molA.setRNABondMobility(BondMobility::Torsion,23-17,25-17);
    molA.setRNABondMobility(BondMobility::Rigid,26-17,29-17);
    molA.setRNABondMobility(BondMobility::Torsion,30-17,35-17);
    molA.setRNABondMobility(BondMobility::Rigid,36-17,39-17);
    molA.setRNABondMobility(BondMobility::Rigid,40-17,45-17);
    cout<<"about to adopt molA"<<endl;
    system.adoptCompound(molA);

	
    GeneralForceSubsystem forces(system); //this is needed only for the VanderWallSphere in WaterDroplet
    HuntCrossleyContact myHuntCrossleyContact(system);
    Force::Custom(forces, new CentrifugalForce(matter));
    system.modelCompounds();        
    State state = system.realizeTopology();
    system.realize(state,Stage::Position);

    

    //DuplexMaker(molA,17-17,22-17,40-17,45-17 ,forces, matter);
    //DuplexMaker(molA,26-17,29-17,36-17,39-17 ,forces, matter);
    //DuplexMaker(molC, molD,forces, matter);
// start adding spheres

    //AddHuntCrossleySpheres(molA,myHuntCrossleyContact, matter);
// done with spheres

 cout<<"check 2.0"<<endl;
     if (! myRigid)
     {
    	//const SimTK::Transform& mytransform2 = matter.getMobilizedBody(molA.getAtomMobilizedBodyIndex(Compound::AtomIndex(0))).findBodyTransformInAnotherBody(state,matter.getMobilizedBody(molA.getAtomMobilizedBodyIndex(Compound::AtomIndex(molA.getNAtoms()-1))));
    	cout<<"check 2.5"<<endl;
    	//Constraint::Weld myConstraint(matter.updMobilizedBody(molA.getAtomMobilizedBodyIndex(Compound::AtomIndex(0))),Transform(Vec3(0)) , matter.updMobilizedBody(molA.getAtomMobilizedBodyIndex(Compound::AtomIndex(molA.getNAtoms()-1))),mytransform2      );
    }
    cout<<"check 3.0"<<endl;

#ifdef CREATE_VTK_WINDOW
    VTKVisualizer display(system, 0.1);
#endif

    VelocityRescalingThermostat * myVelocityRescalingThermostat = new VelocityRescalingThermostat(system,  myTemperature, myReportingInterval);
    ofstream output("mymovie.4.pdb");
    PeriodicPdbWriter * myPeriodicPdbWriter = new PeriodicPdbWriter(system,output,myReportingInterval);
    cout<<"check 4.0"<<endl;
    system.updDefaultSubsystem().addEventHandler(myVelocityRescalingThermostat);
    system.updDefaultSubsystem().addEventReporter(myPeriodicPdbWriter);
    cout<<"check 5.0"<<endl;


    dumm.setCoulombGlobalScaleFactor(myCoulombGlobalScaleFactor);	
    dumm.setBondStretchGlobalScaleFactor(myBondStretchGlobalScaleFactor) ;
    dumm.setGbsaGlobalScaleFactor(myGbsaGlobalScaleFactor);	
    dumm.setVdwGlobalScaleFactor(myVdwGlobalScaleFactor);	
    dumm.setUseMultithreadedComputation(myUseMultithreadedComputation);
    Force::Custom(forces, new CentrifugalForce(matter));

    RungeKuttaMersonIntegrator study(system);
    cout<<"check 6.0"<<endl;


    //scf put back in later
    system.realizeTopology();

    int p = 0;
    stringstream ss1;
    ss1<<"out."<<p<<".pdb";
    ofstream outputFrame(ss1.str().c_str());

    cout<<"p = "<<p<<endl;
    for (int i = 0; i< system.getNumCompounds(); i++)
		{
		cout<<"printing out Compound: "<<i<<endl;	
		//system.updCompound(Compound::Index(i)).writePdb(study.getState(),outputFrame,Vec3(0));
		}
    cout<<"done printing out minimized structure"<<endl;	
    cout<<"multithreading set to : "<<dumm.getUseMultithreadedComputation()<<endl;
    cout<<"Gbsa set to : "<<myGbsaGlobalScaleFactor      <<endl;
#ifdef CREATE_VTK_WINDOW
    display.report(study.getState());
#endif

    TimeStepper ts(system,study);
    ts.initialize(system.getDefaultState());
    ts.stepTo(maxTime);
    return(0);
}

